CHEBYSHEV EXPANSIONS FOR 
SOLUTIONS OF LINEAR DIFFERENTIAL EQUATIONS 



ALEXANDRE BENOIT AND BRUNO SALVY 



Abstract. A Chcbyshcv expansion is a scries in the basis of Chebyshev poly- 
nomials of the first kind. When such a series solves a linear differential equa- 
tion, its coefficients satisfy a linear recurrence equation. We interpret this 
equation as the numerator of a fraction of linear recurrence operators. This 
interpretation lets us give a simple view of previous algorithms, analyze their 
complexity, and design a faster one for large orders. 



1. Introduction 
Chebyshev series are series of the form 

oo 

(1) f{x)^^+Y,CnTn{x), 

^ n=\ 

where T„ denotes the nth Chebyshev polynomial of the first kind. These polyno- 
mials can be defined by 

(2) Tn{cose) ^ cos{ne), 

so that these series behave like Fourier series. Thus in particular, this series con- 
verges pointwise to / on [—1, 1] if / is continuous there, while the convergence is 
uniform if / satisfies a Dini-Lipschitz condition or is of bounded variation (and a 
fortiori if it is differentiable), see, e.g. [8, 11]. Then truncations of the series pro- 
vide polynomials with good approximation properties on the interval [—1, 1], which 
makes these series an interesting data structure for real functions [16]. 

Orthogonality of the T„ leads to the following integral representation of the 
coefficients: 

TT J„l VI - 

We say that / admits a Chebyshev expansion c„T„ when these integrals con- 
verge, the symbol ^ accounting for the factor 1/2 in front of cq in (1). 

In the frequent case when / is a solution to a linear differential equation, Clcn- 
shaw [5] has given a numerical scheme to compute the coefficients without com- 
puting all these integrals. In that case, the coefficients c„ obey a linear recurrence 
equation. A method for the computation of this recurrence has been showed by sev- 
eral authors, first for small orders [6, 10], then in more generality by Paszkowski [13] 
and in the context of (early) symbolic computation by Geddes [7]. We call this 



Key words and phrases. Chebyshev series. Ore polynomials. 

This work was supported in part by the Inria-Microsoft Research Joint Centre. 

To appear in the proceedings of ISSAC'09. 



1 



2 



ALEXANDRE BENOIT AND BRUNO SALVY 



method "Paszkowski's algorithm" . The use of this recurrence to compute the co- 
efficients numerically is discussed in [19]. Paszkowski's method has been further 
improved by Lcwanowicz [9] who gave an algorithm computing a smaller order re- 
currence in some cases. However, Lcwanowicz's algorithm is not much discussed in 
the literature since it looks complicated (see the original article and the comment 
in [19, p. 186]). More recently, other methods have been given by Rebillard [15] 
and Rebillard and Zakrajsek [14]. 

In this work, we give a simple unified presentation of most of these algorithms, 
and design a faster one for large orders. Postponing the proofs and rigorous def- 
initions, the basic idea can be presented by analogy with the computation of a 
recurrence for coefficients of Taylor series. The monomial basis A/„(x) = a;" satis- 
fies 

(3) xMnix) = Mn+iix), M;(x) = nM„_i(x). 

The analogous relations on the Chebyshev polynomials are easily derived from (2) 
and trigonometry: 

(4) 2xT^{x) = Tn+i{x) + r„_i(x), 

(5) 2(1 - x^)T'^{x) = -nTr,+i{x) + nTn-i{x). 

Given a series f{x) ~ ^c„A/„(a;), (3) leads to expressions for the coefficient of 
Mn{x) in xf and /': multiplication by x maps to a negative shift on the indices; 
differentiation maps to a positive shift of the index followed by multiplication by n+ 
1. Algebraically, we thus get an algebra morphism mapping x to X and d/dx to D, 
with X := S~^, D := (n-|-l)5'. Here, S denotes the shift operator: u{n) h^- u{n+l), 
that does not commute with multiplication by n. Now, if / is solution of a linear 
differential equation 

Pkix)f^''Hx) + ---+po{x)fix)^0, 
we deduce a recurrence operator pk{X)D^ + ■ • ■ + pq{X) for its Taylor coefficients. 

Example 1. The simplest example is the exponential, for which f — f = trans- 
lates into D — 1 = (n + l)S'— 1 (1 denotes identity), which gives the recurrence 
(n -|- l)c„+i — c„ = satisfied by c„ = 1/n!. 

The procedure for a series f{x) = J^'^nTnix) starts similarly: multiplication 
by X maps to 

(6) X ■.= {S + S-^)/2. 

The difference comes from the factor (1 — x^) in (5). The operation of differentiation 
followed by multiplication by (1 — x^) is readily seen to map to {S — S~^)n/2, but no 
simple linear operation for the Chebyshev coefficients of /' exists. The idea at this 
stage is to divide by 1 — a;^ afterwards by introducing a formal inverse of 1 — X^. 
Thus we write D :— (1 — X'^)^^{S — S^^)n/2. This can be further simplified since 
1-X'^ = -{S- 5-1)2/4, so that 

(7) D ■.= 2{S-^ - S)-^n. 

We call such an expression a fraction of recurrence operators. 
Example 2. For the exponential, wc now get 

D-l = 2{S-^ - Sy^n - 1 ^ {S-^ - S)-\2n - (5"^ - S)). 
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The last term is an analogue of reduction to the same denominator. The final factor 
will be called the numerator of the fraction. It corresponds to the recurrence 

(8) 2nc„ - c„_i + c„+i = 0. 

It turns out that in this example, the Chebyshev coefficients are known: c„ = 
2J„(1), where J„ is the modified Bessel function of the first kind, and they do 
satisfy (8). 

This example generalizes. We show here that all the algorithms mentioned above 
can be interpreted as first rewriting the input linear differential equation in one way 
or another, then applying the morphism above and finally returning the numerator 
of the result. In the case of Lewanowicz's algorithm the fraction is normalized (its 
numerator and denominator are relatively prime), which is why its output may 
have smaller order. 

In Section 2, we give the formal setting for fractions of recurrence operators, to- 
gether with the basic algorithms. This is then applied to the specific case of Cheby- 
shev series in Section 3. Then we give a compact presentation of Paszkowski's and 
Rebillard's algorithms, provide a complexity analysis and design a faster algorithm 
in Section 4. We briefly comment on the different approach taken by Rebillard and 
Zakrajsek in §4.5. We conclude in Section 5 with a few examples. 

2. Fractions of Recurrence Operators 

We use Ore's framework of non-commutative polynomials [12], that we now 
recall. 

2.1. Ore Polynomials. The rings of linear differential operators and of linear 
recurrence operators are special cases of rings of Ore polynomials. They possess 
the commutation rules 

-^p(x) = p{x)^ + p'{x); Sp{n) = p{n + 1)5. 
ax ax 

More generally, a ring of polynomials in an indeterminate d with coefficients in a 
field K is an Ore polynomial ring when its product is defined by associativity from 

(9) dp^a{p)d + S{p), peK 
where for all a and b in K, 

(T{a + b) ^ (T{a) + a{b), a{ab) ~ (T{a)a{b), a{a^^) ^ a{a)~'^ , 
6ia + b) = 6{a) + S{b), S{ab) = a{a)d{b) + S{a)b. 

The ring is denoted K{d;a,S). Linear differential operators are obtained with a = 
Id and S = d/dx; linear recurrences operators with a = S and i5 = 0. 

The main property of these rings is that the degree (with respect to d) of a 
product is the sum of the degrees of its factors. (In particular, there are no zero- 
divisors). From there, it is not difficult to write an algorithm for Euclidean division 
on the right. Once right Euclidean division is available, the Euclidean algorithm 
and its extended version follow and can be used to compute: greatest common right 
divisors, denoted gcrd; least common left multiples, denoted Iclm; the corresponding 
cofactors for the Bezout identity and for the Iclm [12, 4]. 

When a is invertible, we also get Euclidean division on the left, and from there 
greatest common left divisors (geld), least common right multiples (Icrm) and the 
corresponding cofactors by the Euclidean algorithm. If moreover 6 = 0, as is 
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the case for recurrence operators, it is also possible to define Laurent polynomials 
with 99^^ = d^^d = 1 and d^^a — a^^ {a)d^^ . These are denoted K{d,d~^;a). 

The rings wc use in this work arc the ring of linear differential operators denoted 
Q{x){dx; Id, d/dx) and the ring of linear recurrence operators Q{n){S, S^^; S) (with 
a different meaning for both S). 

Apart from their non-commutativity. Ore polynomials generally behave like or- 
dinary polynomials. A notable exception is divisibility. 

Example 3. The recurrence operator P — [n + 1)^^{S + I) is relatively prime 
with Q = nS + n + 2, but Q is a right divisor of P^. 

Still, the following property holds (and similarly for gcrd's when they exist): 

(10) gc\d{AB,AC) = Agc\d(B,C). 

Indeed, A is a left divisor of gcld(^i?, AC) and the remaining factor has to be a 
left divisor of both B and C. The converse divisibility is clear. 

In order to distinguish the action of an operator from the product in these rings 
of operators, which corresponds to composition of actions, we use the notation • for 
the former. Thus dx ■ f = f , S ■ Un = Un+i. 

2.2. Fractions. Ore's construction of fractions parallels the commutative case. 
Given two non-zero polynomials Qi and Q2, by definition of the Iclm, there ex- 
ist two polynomials Qi and Q2 such that 

lclm(Qi,02) = QiQi =Q2Q2. 

With this notation, the pairs {Pi,Qi) and (^2,^2) arc called equivalent when 
QiPi = Q2^2- This can be verified to be an equivalence relation and the class is 
called a fraction and denoted Qi^Pi (which is equal to Q2^P2)- This construction 
makes the set of fractions a non-commutative field. 

Reduction to the same denominator for sums is given by 

(11) Qi'Pi + Q2^P2 = \c\in{Qi,Q2)-\QiPi + Q2P2), 

as can be checked by left multiplication with lclm((3i, Q2)- 

To compute the reduction of a product of two fractions Q^^Pi, Q2^P2, the start- 
ing point is the Iclm of Q2 and the numerator Pi. There exist two polynomials Pi 
and P2 such that 

\c\m{Q2,Pl)=PlPl=Q2Q2. 

Then, Q^^Pi = [PiQiY^PiPi and Q^^P2 = (Q2Q2)"'Q2P2, so that finally 

(12) Qi^PiQ2^P2 = {PiQiy^Q2P2- 

2.3. Irreducible Fractions. Having in mind our use of fractions for recurrence 
operators, wc now concentrate on the case when a is invcrtible, so that gcld's are 
available. The results here are probably known, but we did not find them in the 
literature. 

A fraction Q^^P is called irreducible when gcld(P, Q) = 1. 

Proposition 1. Assume a is invertible and let A~^B be a fraction. Then there 
exists an irreducible fraction equal to A~^B. Moreover, its numerator and denom- 
inator are unique up to a factor in K. 
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Proof. Existence follows from dividing out numerator and denominator by gcld(A, B). 
Assume Qi^Pi = Q2^P2 and gcld(Pi,(5i) = gcld(P2,Q2) = 1- By definition of 
equivalence, QiPi = Q2P2, where QiQi = Q2Q2 ~ lclm(Qi, (52)- Moreover this 
Iclm relation implies gcld(Qi, Q2) = 1- Now, 

Qi = Qi gcld(Pi, Qi) = gcld(giPi, giQi) 

= gcld(g2P2, Q2g2) = Q2 gcld(P2, Q2) = Q2, 

where we use (10). But since gcld(Qi, Q2) = 1, necessarily Qi = Q2 = 1 and then 
Pi = P2 and Qi = Q2. □ 

The following lemma is useful in the computation of recurrences for Chebyshcv 
series. 

Lemma 1. Assume a is invertible and let Q2^P2 be an irreducible fraction and Pi 
a polynomial. Then P^^Q2P2 with Q2Q2 = Pi Pi = \c\m{Q2, Pi) is irreducible and 
equal to P1Q2 Pj- 

Proof. We have gc\d{Q2, Pi) = 1 by definition of the Iclm. The polynomial g = 
gcld(Pi, g2P2) is a left divisor of Pi Pi = Q2Q2, and therefore is a left divisor of 

gcid(g2P2, g2g2) - g2 gcid(P2, g2) = g2. 

Thus g is a left divisor of both Pi and g2, hence is 1. □ 



3. Recurrences for Chebyshev Coefficients 

We now have the theoretical tools to prove that a morphism from linear differen- 
tial operators to linear recurrence operators produces fractions whose numerators 
give recurrences for the coefficients of Chebyshev series solutions. 

The algorithms then become easy to state, their algorithmic difficulty being 
concentrated in the Euclidean algorithm in the previous section. 

3.1. Morphism. We define a morphism of Q-algebras from Q[x]{dx','iii, d/dx) to 
the field of fractions of Q{n){S, S~^; S) by 

fix) =X:=^(S + S-'), f{dx) = D := [S-^ - S)-\2n). 

The proof that is a well-defined morphism of non-commutative rings reduces to 
checking the commutation ip(dxx) = ip{xdx + 1). Indeed, 

XD + l = ]^{S + S-^){S-^ - S)-^{2n) + 1 
= {S-^ - S)-^{S + S-^)n + \ 

= (5-1 - S)-^({{n + 1)S + (n - l)S-^) + (5-1 - 5)) 
= DX. 
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Algorithm 1 Lcwanowicz' algorithm 

Input: L := Ej=oP»(^)*5' 

Output: (P,Q) such that Lp{L) = Q^^P 

Q := 1 

for all i from /c — 1 to do 

Compute lclm((5-i - S), P) ^ PP = U{S-^ - S). 
Q:=PQ 

P := U2n + Qp,{X) 
end for 
return (P, Q) 



3.2. Horner's Rule and Lewanowicz' Algorithm. 

Proposition 2. Let L — pi^{x)d^ + • • • + dxPo{x) he a linear differential operator 
in <Q\x\{dx\^A,d/dx) . The evaluation of ip{L) by Horner's rule 

ifiL) = {■ ■ ■ {pk{X)D + pk-i{X))D + ■ ■ ■)D + po{X) 

using Eqs. (11) and (12) for the computation of sums and products produces a 
fraction Q^^P that is irreducible. 

The algorithm deduced from this statement (Algorithm 1) is due to Lewanowicz. 
It is made very clear by the use of fractions of recurrence operators. The proof that 
the numerator of its output gives a recurrence for the Chebyshev coefficients is 
given in the next section. 

Proof. We prove that each iteration of the loop produces (P, Q) that are relatively 
prime and such that 

(13) Q-'P = ^(Pfe(x)9^-* + • • • + p,). 

Initially, i = k and f{pk{x)) is a polynomial, so that Q = 1 and the property holds. 
If it holds for Mi, the next stage of the loop computes Q^^PD + pi-i{X). Recall 
that D = (5-1 - 5)-i(2n). Then let lclm((5-i - S),P) = PP = U{S-'^ - S). 
Lemma 1 applied to the inverse {S^^ — S)P^^Q implies that gc\d{PQ,U) = 1. It 
follows that gcld{PQ,U + PQpi^i{X) / {2n)) = 1. Again by Lemma 1 apphed to 
the inverse, multiplying by 2n on the right preserves irreducibility and the property 
holds for Af,_i. □ 

Wc quote without proof the following result. 

Proposition 3 (Lewanowicz). When the leading coefficient pk{x) of the differential 
equation does not vanish at 1 or —1, then all the gcrd's are trivial, Q = Z)~' at 
step i and the resulting Q is D^^ . 

This is related to the fact that 4(^2 - 1) = [S-^ - Sf. 



3.3. Chebyshev Expansions. 
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3.3.1. Main Theorem. We now prove our main result: the morphism defined above 
behaves as expected with respect to Chebyshev expansions. 

Theorem 1. Let L = ^0(2;) + ■•• + Pk{x)d^ be a linear differential operator of 
order k with polynomial coefficients. Let / G C'''(] — 1, 1[) be such that either of the 
following hypotheses holds: 



(H) / —j^=^= dx is convergent; 

J — 1 V J- X 



1-1 VT 

(1 - x^)'' f'-'^Hx) 2 ■ 

(H') J ^ dx is convergent and (1 — x y\pi, i = 0, . . . , k. 

Then f admits a Chebyshev expansion UnTn, L ■ f admits a Chebyshev expan- 
sion ^* w„r„ and the sequences (w„) and (w„) are related by P ■ Un = Q ■ w„, for 
any {P, Q) such that Q^^P = ip{L). In particular, if L ■ f ~ Q, then the Chebyshev 
coefficients of f satisfy P ■ Un = for any numerator of ip{L). 

The easy case is when (H) holds. Hypothesis (H') makes it possible to deal with 
some functions that are singular at ±1, but whose singularities are not "too bad": 
they are regular singular points. 

Proof. First, convergence of the integral in (H) or (H') implies convergence of the 
analogous integral where /'^'^^ is replaced by for i — fc — 1, . . . , as well as the 
integrals where these functions are multiplied by r„(x), n £ N. This shows that 
both / and L ■ f admit Chebyshev expansions. 

If the result holds for any numerator of ip{L) then in particular it has to hold 
for the numerator of its irreducible form. Conversely, if P • m„ = Q ■ Vn, then 
RP ■ u„ = RQ ■ Vn for any R, so that it is also sufficient to prove the result for an 
irreducible 4'{L). 

Lemma 2 (Basic Cases). Under the same hypotheses, the result holds for L a 
constant times identity, L = x, L = dx if (H) holds, L = (1 — x^)dx if (H') holds. 

Proof. If i = A is a constant times identity then P — X, Q = 1 and u„ = Au„ 
clearly holds. 

If L = X, Eq. (4) implies 

2 f' f{x)xTn{x) 2 /(x)(r„+i(a;)+r„_i(x)) ^ 

■ dx ~ — , dx — X ■ Un. 



li L — dx and (H) holds, we use the following variant of Eq. (5) when n 7^ 



V 2nVl - a;2 / n^^ 



VT- 

that can be checked from (2). The continuity of /' and the convergence of the 
integral in (H) imply that integrating by parts is possible and this gives 

2 f{x)Tn{x) 



■ dx 



2/(x)yr^T^ 



2 f'{x){Tn-l{x)~Tn+l{x)) 



-1 



dx^D 



-1 



Both limits of the term between brackets are 0, by convergence of the integral u„. 
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The case when n = reduces to checking v-i = vi, that does not depend on /. 
If L = (1 — x^)dx and (H') holds, we start from 

An argument similar to the previous one then gives 

(n + l)u„+i - (n - l)u„_i = - / /(a:) , dx 

= 2- (1 - x')fix)^I^dx ^ 2t;„, 
which proves the result since ^((l - x'^)dx) = (1 - X'^)D = (5 - 5"i)ri/2. □ 



Lemma 3 (Product). Assume the result holds for L2 with f , as well as for another 
operator Li with L2 ■ /• Let ip{Li) = Qi^Pi and 1^9(^2) = Q2^P2, these fractions 
being irreducible. Let lclm(Q2, ^'i) = F'lF'i = Q2Q2, and assume (PiQi)~^Q2P2 is 
irreducible. Then the result holds for L1L2 with f . 

Proof. Let w„, Wn, m„ be related by Pi • u„ = Qi • w„, P2 ■ Un ^ Q2 ■ Vn. Then 

PlQl ■ Wn = A-Pl • Vn = Q2Q2 ■ Vn Q2P2 ' M„, 

whence the result. □ 



As a consequence, the result holds when L = Ax' is a monomial, by induction. 

Lemma 4 (Sum). Assume the result holds for an operator L with f and for a 
polynomial p with the same f. Then it holds for L + p with f. 

Proof. Let (piL) — Q^^P be irreducible. If P • u„ = Q • z;„, Wn — p{X) ■ u„, then 

Q-{Vn+ Wn) = (P + Qp{X)) ■ Un. 

This proves the property for L +p since gcld(Q, P + Qp{X)) = gcld((5, P) = 1. □ 



The result now holds for L an arbitrary polynomial, as a sum of its monomials. 

Let finally 6' = 9^ if (H) holds and 6* = (1 — x^)dx if (H') docs. In both cases, 
L can be written qk{x)9^ + ■ ■ ■ + qQ{x) with polynomial q^, . . ■ , qo. The hypothesis 
on / implies that the result holds for L = 1 with 6'' • / for i = 0, . . . , /c and therefore 
also for L = qi{x) with 6'' • / by Lemma 3. 

Let Lk = qk{x) and Li = Li+iO + q,; for i = fc — 1, . . . , 0. Let 0(ii) = Q^^Pi 
be irreducible. We prove by induction that the result holds for Li with 9^f. For 
i = k, the result has just been proved. If the result holds for P^+i with 9^^^f, 
then we obtain an irreducible (5~|_\Pi+i0: when 9 = dx this follows from Lemma 3, 
while when 9 = {\ ~ x'^)dx, ij+i itself is a polynomial (by induction). Thus the 
result holds for Li+i9 with 9^ f . Since it also holds for qi{x) with 9'' f and qi{x) is 
a polynomial, wc get the result for their sum by Lemma 4. Thus by induction the 
result holds for Lq with /, which concludes the proof of Theorem 1. □ 
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3.3.2. Examples. 

Example 4. The function cxp(a;) satisfies (H) for any k. This proves the recur- 
rence (8) computed in Example 2. 

Example 5. The function (1 — x^)~^^^ is annihilated by 2(1 — x^)dx — x. Hy- 
pothesis (H) does not hold, but (H') does. Application of the morphism gives 
P = {2n + 3)S'2 - (2?! + 1), Q = -25", so that the theorem asserts that the Cheby- 
shev coefficients satisfy 

(14) (2n-H3)c„+2 = (2n+l)c„. 

The actual values can be computed by standard properties of the Beta integrals 
and indeed 

if ri is odd, 



Cri 



V^r(t+|) otherwise. 



Example 6. The function arccosx gives an example showing that analytic hy- 
pothesis such as (H) or (H') are necessary. This function is annihilated by L = 
(1 — x^)d^ — xdx- Direct application of the morphism gives P = n^, (5 = 1, which 
would suggest that the recurrence is n^Cn = 0. However, neither (H) nor (H') holds 
in this case. Left multiplying L by (1 — a;^) gives a new operator such that (H') 
holds. Then the theorem proves that the coefficients are annihilated by 

(15) {n + 4)^S^ -2{n + 2f + n^. 
This can be checked against the actual coefficients: 

if 71 = 0, 

if n > is even, 
otherwise. 

4. Algorithms 

We now cast the algorithms of Paszkowski [13] and Rebillard [15] as computations 
of the numerator of a fraction of recurrence operators. We also propose a new faster 
algorithm. All three algorithms compute the same recurrence. Starting from 

k 

(16) L = ^p,(a;)a;, 

1=0 

they avoid the need for fractions by replacing differentiations by integrations, ex- 
ploiting the polynomial 




I -D-^ = { — ] (-S + S 



These algorithms compute the polynomial l''(p{L), that is a numerator of (p{L) = 
I~^I*'tp[L). Thus, by Theorem 1, their result is a recurrence operator annihilating 
the coefficients of Chcbyshcv series solutions of L. 

If pa;(1)pa;(— 1) 7^ 0, Proposition 3 shows that is the denominator of the 
irreducible fraction and therefore in that case all algorithms compute the irreducible 
fraction. Otherwise, the result of these algorithms may have larger order than that 
returned by Lewanowicz' algorithm. 
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Example 7. The function (1— x^)^^/^ has been deah with in Example 5. Lewanow- 
icz' algorithm returns the second order recurrence (14). The numerator returned 
by the other algorithms has order 4: 

(2n + l)c„ - 4(n + 2)c„+2 + (2n + 7)c„+4 = 0. 

It is however possible to recover the smaller order recurrence: the geld oi A = 
(2n + 7)S'4 - 4(n + 2)S^ + {2n + 1) with / is /, so that A factors as A = (5^ - 1)P 
with P as in Example 5. 

More generally, dividing out the result of the computation of I^(f){L) on the left 
by the geld with yields the result of Lewanowicz's algorithm. 

4.1. Paszkowski's Algorithm. The starting point of Paszkowski's algorithm is 
to rewrite L from (16) as 

k 
1=0 

The polynomials qi can be computed inductively starting with qt — Pk and sub- 
tracting dl^qk to produce a smaller order operator. Then 

k 

(17) /V(i) = E^'"^9.(^)- 

i=0 

Algorithm 2 follows. 



Algorithm 2 Paszkowski's Algorithm 

Input: L :=^J2i=oP^i^)^x 
Output: l''ip{L) 

Compute g-i's such that L ~ X]i=o (^x1i{^) 

R:=qk{X) 

for all i from 1 to fc do 

R:= R + Pqk-^{X) 
end for 
return R 



4.2. Rebillard's Algorithm. The starting point of Rebillard's algorithm is the 
identity 

Xk = I'^XD'' = (2n)-i((n + k)S + {n - k)S-^), 
that follows from an easy induction. From there, he deduces 

k k 

/V(i) = 5^/V.(^)^'/'"'' - ^p.(x,)/'=-\ 

i=0 i=0 

Algorithm 3 follows. 
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Algorithm 3 Rcbillard's Algorithm 

Input: L = Y!l=oP^{x)dl 
Output: I^ifiL) 

Compute pi{Xk), i = . . . ,k 

R:=Pk{Xk) 

for all i from 1 to k do 

R■.= R+Pk-^{Xk)P 

end for 
return R 



4.3. Complexity Analysis. Wc now give a complexity analysis of Paszkowski's, 
Rcbillard's and Lcwanowicz' algorithms. This reveals a source of inefficiency for 
large orders, that we correct in our new algorithm in the next section. 

We need to consider the size of polynomials in two variables n and S. We say 
that a polynomial has bidegree {m,p) in {n, S) when it has degree m in n and p 
in S. 

First, we state more precisely the shape of /*. 
Proposition 4 (Rebillard [15]). For all i e N*, 

1 




k=l 

where r{i) = 2'n'[llZ,\{n^ - P), 



s(fc) = (-1)'= in-i + 2k){n + k+ - z + l)fe„i, 

and we use the Pochhammer symbol {a)i = a{a + 1) • • • (a + i — 1). 

In particular, the bidegree of r{i)P in (n, S) is (i — 1, 2i). The proof is a tedious 
but easy induction that we omit here. From this formula follows a precise estimate 
of the size of the polynomial we are computing. 

Corollary 1. If L in (16) has bidegree {d,k) in {x,dx), then r{k)I^Lp{L) is a 
polynomial of bidegree in (n, S) at most {2k ~ 1, 2(fc + d)). 

Proof. First, L can be rewritten as in Paszkowski's algorithm X]i=o with 
deg Qi < d. The identity 

shows that this is a polynomial in n. Each term of the sum is the product of a 
polynomial of bidegree (2(fc— z),0), apolynomial of bidegree (i— l,2i), apolynomial 
of bidegree at most {0,2d). Thus each summand has bidegree at most {2k — i — 
l,2i + 2d), whence the result. □ 



Proposition 5. Given L as above for input, Paszkowski 's algorithm requires 0{dk^) 
arithmetic operations in Q. 
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Proof. The first step is the computation of the g^'s from the pi's. The inductive 
method requires only 0{dk^) arithmetic operations. Using ideas from [3], it is 
actually possible to decrease this complexity further to 0{M{dk)) operations [2] 
(here, M is the complexity of polynomial product, see, e.g.. [18]). 

The next step is the loop. The main cost in step i is the multiplication of /* by 
qk-i{X). We multiply a polynomial of bidcgrce {i — 1, 2i), with a polynomial in S 
only, of degree 2d. The cost of this multiplication is 0{i^d) arithmetic operations. 
Summing for i up to fc gives the result. □ 

Proposition 6. In the same conditions, Rebillard's algorithm requires 0{d^k + 
d'^k^) arithmetic operations. 

Proof. The first step is the computation of the Pi{Xk). The polynomial has 
bidegree (3i,2i) in {n,S). Then ea.ch. pi{Xk) can be computed in 0{d^) operations 
and all of them in 0{d^k) operations. 

The cost of the zth step of the loop is dominated by the cost of the multiplication 
oi pk-i{Xk) by P. The polynomial has bidegree (3o?, 2c?) in {n,S), while 

P has bidegree {2i — l,2i). Naive multiplication then requires 0{d^i^) operations. 
Summing over k gives the result. □ 

The output of Lewanowicz' algorithm is different in general. We give a compar- 
ison in the cases when it coincides. 

Proposition 7. In the same conditions, and if all the gcrd during its execution are 
trivial, Lewanowicz' algorithm requires 0{dk^) arithmetic operations. 

Proof. We only give a sketch. When all gcrd's are trivial, it turns out that the 
computation of Iclm's and cofactors is of the same order of complexity as the com- 
putation of the product Qpi^ where moreover Q = This is the same as in the 
analysis of Paszkowski's algorithm. □ 

4.4. New Fast Algorithm. Wc now give another algorithm for the same operator 
l'^(p{L). The design of our algorithm is motivated by computational complexity is- 
sues. In the analyses above, most of the complexity comes from the fact that during 
the computations, the bidegrees of the intermediate polynomials grow linearly and 
they are multiplied by polynomials of fixed degree. Instead, we aim at balancing 
degrees so as to make use of the recent fast algorithm for the product of linear 
differential operators [17, 3], that we denote FFT-mult. We achieve the following 
complexity. 

Theorem 2. Algorithm 4 computes the recurrence operator I^Lp{L) in 0{{d + 
k)k'^^^) arithmetic operations. 

Here, w is a feasible exponent for matrix multiplication with coefficients in Q 
(see, e.g., [18]). We now prove this result. 
Starting from (17), we write 

k (-1 k 

J2 I\k-^{X) = + E I'"'ik-^{X) ~ P(o,.. + /'Po,.. 

i=0 i=0 i=t 

We choose I = \k/2'\ and apply the same idea recursively. 

The time consuming part of the computation is the product I^P(e....,k)7 for which 
we give a specialized algorithm. 
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Algorithm 4 Divide and Conquer 

Input: Polynomials ao(a;), . . . , ak{x) 

Output: P(o,...fc) =Eto-^'«»(^) 
if fc = then 

return ao(X) 
else 

£ := [fc/2] 

Compute recursively P[o,...j-i) and P(e,...,k)- 
return P(o,.../-i) + I^P{i,...,k)- 
end if 



Algorithm 5 Fast Multiplication 

Input: P and P := ELo-^'a»(^) 
Output: /''P 

Decompose P as in Eq. (18) 

R:=0 

for all i from to [{£ + d)/£\ do 

R:=R + FFT-mult(r(^)/^ 
end for 

return l/r{2l)R. 



To simplify the presentation, assume /c = 2i. Corollary 1 implies that has 
degree 2£ in S, P{i^...^k) has degree at most 2£ + 2d in S*. They have rational function 
coefficients whose degrees are also bounded by this result. If d is large, the degrees 
in S arc unbalanced, so we first decompose 

(18) P(,,...,fc) = Ao{n, 8)8-"-' + Ai(n, 8)8~''+'+^ + • • • , 

where the AiS have degree at most 21 in 8. Note that this decomposition is only 
an extraction of coefficients and does not use any arithmetic operation. 

We are thus left with the product of with the A^'s. Although both have ra- 
tional function coefficients and thus cannot be multiplied directly by FFT-mult, we 
also have that r(£)J^P(£_...^fc) has polynomial coefficients in n of degree at most fc— 1 
and therefore so does r{t}I^Ai. To perform the product efficiently, we make use 
of the fact that FFT-mult proceeds by evaluation and interpolation: during the 
evaluation phase, we evaluate the rational function coefficients of Ai as if they 
were polynomials (and within the same complexity thanks to our degree bounds), 
avoiding the zeros —I, . . . , £ of their denominators; similarly, we evaluate the poly- 
nomial coefficients of r{£)I^. Then we compute the necessary products. With the 
bounds on the degree in n we have for the polynomial coefficients in the result, 
the interpolation phase then returns the result. The complexity of each of these 
multiplications is thus 0(£") operations. The algorithm for this multiplication is 
summarized in Algorithm 5. Note also that a constant factor can be saved by not 
recomputing the "FFT" of r{£)I^ at each time. 

Proposition 8. The cost of multiplying by X]i=o P'^ii-^) with dega^ < d using 
Algorithm 5 is 0{{£ + d)£'^~^) arithmetic operations. 
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Proof. We have seen that each muhiphcation r{i)I^Ai has complexity 0{i'^). This 
is performed {£ + d)/£ times. Right muhiphcation by powers of S does not use 
any arithmetic operations. The additions require a smaher number of operations, 
whence the resuh. □ 

Now, let r(fc, d) be the complexity of Algorithm 4. Using this proposition, we 
get 

T{k, d) = 2T(fc/2, d) + 0{{d + fc)fc"-i), 

the complexity estimate in Theorem 2 follows from the convergence of the geometric 
series. (Again, a constant factor can be saved by computing the powers of / only 
once.) 

4.5. Algorithm by Rebillard and Zakrajsek. In [14] . an algorithm of a different 
nature is proposed. It docs not compute a numerator of (t>{L), but manages in 
some cases to derive a smaller order recurrence corresponding to a right factor 
of the numerator of 4>{L). We plan to come back to this algorithm in connexion 
to minimality issues in future work. Here, we merely give a few indications and 
comments on special cases. 

Example 8. The following is taken from [14]. Starting from the differential opera- 
tor L — {x+ l)^i9^ — {x+ l)dx + X + 7/4, the computation of by Lewanowicz's 
algorithm leads to a numerator of order 4, whereas the algorithm in [14] produces 
one of order only 3. We note that this operator can also be obtained by Lewanow- 
icz's algorithm, applied to dxL instead of L. In many cases, this technique applies. 

Since the algorithm in [14] computes a right factor of the numerator of 
analytic hypotheses such as (H) or (H') in our Theorem 1 are necessary. 

Example 9. In the case of arccos(x), the numerator of (f>{L) is a constant (see 
Example 6), so that this is also the result of the algorithm in [14]. Starting from 
(1 — x'^)L, we obtained an operator of order 4 in Ex. 6. On this operator, the 
algorithm in [14] returns an operator of order 2. 



5. Examples 

The fast algorithm does not lend itself easily to an efficient implementation in 
Maple, since it relies on fast evaluation/interpolation and fast matrix product. 
We have however implemented the slow algorithms in Maple and show how other 
algorithms from computer algebra can sometimes be applied to the resulting recur- 
rences, so that nice expression for the coefficients can be recovered. We have also 
implemented variants of Horner-like evaluations that seem to perform well, see [1]. 

Example 10 (arctan). Starting from (.t^ + 1)9^ + 2xdx, we get 

ncn + (6n -I- 12)c„+2 + (n + 4)c„+4 0. 

The initial conditions are computed by Maple as cq = C2 = and ci = 2^/2 — 2, 
C3 = (14 — 10\/2)/3. The recurrence can then be solved by Petkovsek's algorithm 
and we get 

arctan(x) = 2 V (1 - ^f+i T^k+iix) ^ 
t^o 2fc + l 
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Example 11 (error function). Starting from d'^ + 2xdx, we get a more complicated 
recurrence: 

(n^ + in)cn + {2n^ + 12n^ + 24ri + 16)c„+2 - {n^ + 5n + 4)c„+4 = 0. 
A closed form is known to be 



but it seems that the algorithms in computer algebra are not strong enough to find 
this automatically, yet. 

Example 12 (arctanh). Starting from L = {x^ — 1)9^ + 2xdx, we get 

[n + 2)c„+2 - ncn = 

by computing the numerator of 4>{L). Although neither of our hypotheses (H) or 
(H') holds here, this result is correct, as can be checked from the expansion 



Again, this suggests that more work on obtaining a recurrence of minimal order is 
necessary. 
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